2b.  OECLASSlFlCATION  /  OOWNG 


4.  PERFORMING  ORGANIZATION  REPORT  NUMB 

< 

Technical  Report  No.  26 


W 


S.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 
ONR 


6s.  NAME  OF  PERFORMING  ORGANIZATION 
Oakland  University 


6b.  OFFICE  SYMBOL  7a.  NAME  OF  MONITORING  ORGANIZATION 

<"  Office  of  S.v.1  Research 


6c  ADDRESS  (Gey,  State,  and  ZIP  Code) 
Department  of  Chemistry 
Rochester,  Michigan  48063 


7b.  ADORESS  (City,  State,  and  ZIP  Cod* ) 
800  N  Quincy  St. 
Arlington,  VA  22217 


Ba.  NAME  OF  FUNDING/ SPONSORING 
ORGANIZATION 


8b.  OFFICE  SYMBOL  I  9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 
(If  applicable)  I 


Office  of  Naval  Research 


8c  ADDRESS  (C/ty,  Sutt.ind  ZIP  Cod*) 
800  N.  Quincy  St. , 
Arlington,  VA  22217 


10  SOURCE  OF  FUNO'NG  NUM8ERS 


PROGRAM 
ELEMENT  NO. 


4133006 


II  TITLE  (Includ*  Security  ClauiAcation) 

Monte  Carlo  Simulation  of  Primitive  Atom* transfer  Reactions  in  Solution 


W  PERSONAL  author(S) 

P.  P.  Schmidt  and  Lesser  Blum 


13a.  TYPE  OF  REPORT 

Technical 


Technical  frov 


16.  SUPPLEMENTARY  NOTATION 

submitted  to  J.  Phy.  Chem. 


13b.  TIME  COVEREO 
FROM _  TO 


14.  DATE  OF  REPORT  (Y*»r.  Month.  Oty)  15.  PAGE  COUNT 

August  5,  1987  47 


COSATI  CODES _ I  18.  SUBJECT  TERMS  (Continue  on  reverie  if  neceuary  and  idtntify  by  block  numbtr) 

FIELO  |  GROUP  I  SUB-GROUP  I 

Atom  transfer  (transition  state)  theory 


19.  ABSTRACT  (Continue  on  reverse  if  necessary  and  idtntify  by  block  numbtr) 


MENT  A 
Approved  fox  public  release 
Distribution  Unlimited 


The  objective  of  this  work  is  develop^ a  computational  algorithm 

which  combines  saddle  point  and  Metropolis/Monte  Carlo  optimization  to 
investigate  reactions  in  solution;  the  reactions  involve  atom  transfer  on  an 
adiabatic  potential  energy  surface.  JZe-aeek  -a  me the d  to  calculate  the  value 
of  the  rate  constant  in  the  form  of  the  simple  Arrhenius  equation  for  the 
jump  rate  v , 


v  -  <w>exp(-(<E*>  -  <E(j>)/kT^, 

in  which  <E*>  is  the  average  energy  of  the  transition  state,  <E<>>  is  the 


20.  DISTRIBUTION /AVAILABILITY  OF  ABSTRACT 
□  UNCLASSIFIEO/UNUMITIO  □  SAME  AS  RPT 

□  OTIC  USERS 

22s  NAME  OF  RESPONSIBLE  INOIVIOUAL 

OO  FORM  1473  r.i  • 


21  ABSTRACT  SECURITY  CLASSIFICATION 


22b  TELEPHONE  (Include  Arts  Code)  22c.  OFFICE  SYMBOL 


91  AP9  h*  E.nt  ’  **  *  “ 


I 


g 


I 


Mi 


i 


mmmmm 


average  energy  of  the  initial  state,  and  <to>  is  the  average  frequency  of 
passage  through  the  transition  state.  Individual  configurations 
in  the  Metropolis  sample  allow  either  for  passage  of  the  reactive  species 
over  the  top  of  the  barrier  or  tunnelling  through  the  barrier.  The  average 
frequency  <ui>  reflects  this  situation.  Because  the  Metropolis  sampling 
method  deals  with  discrete  collections  of  particles,  with  specified  forces 
of  interaction,  the  transfer  frequencies  for  over- the -top  of  the  barrier  and 
tunnel  transfers  can  be  determined  in  terms  of  the  actual  interactions  used 
instead  of  using  non-specific,  model  potential  energy  functions  for  the 
barrier . 


OFFICE  OF  NAVAL  RESEARCH 


Contract  No.  N00014-77-C-0239 
R&T  Code  4133006 

Technical  Report  No.  26 


Monte  Carlo  Simulation  of  Primitive  Atom- transfer 
Reactions  in  Solution 


by 

P.  P.  Schmidt  and  L.  Blum 


Prepared  for  Publication  in  the 
J.  Phys.  Chem. 


Oakland  University 
Department  of  Chemistry 
Rochester,  Michigan  48063 


July  1987 


Reporductlon  in  whole  or  in  part  is  permitted  for 
any  purpose  of  the  United  States  Government 


Approved  for  Public  Release;  Distribution  Unlimited 


Monte  Carlo  Simulation  of  Primitive  Atom- transfer  Reactions 


in  Solution 


P.  P.  Schmidt 


Department  of  Chemistry 
Oakland  University 
Rochester,  Michigan  48063 


Lesser  Blum 


Department  of  Physics 
University  of  Puerto  Rico 
Rio  Piedras,  Puerto  Rico 


& 

IS 

I 

M 

i 

I 


$ 

I 

K 

S 

'A 

A 

i 

i 


Abstract 


The  objective  of  this  work  is  to  develop  a  computational  algorithm 
which  combines  saddle  point  and  Metropolis/Monte  Carlo  optimization  to 
investigate  reactions  in  solution;  the  reactions  involve  atom  transfer  on  an 
adiabatic  potential  energy  surface.  We  seek  a  method  to  calculate  the  value 
of  the  rate  constant  in  the  form  of  the  simple  Arrhenius  equation  for  the 
jump  rate  v, 

v  -  <w>exp[-(<E*>-<E°>)/kT]  , 

in  which  <E*>  is  the  average  energy  of  the  transition  state,  <E°>  is  the 
average  energy  of  the  initial  state,  and  <w>  is  the  average  frequency  of 
passage  through  the  transition  state.  Individual  configurations 
in  the  Metropolis  sample  allow  either  for  passage  of  the  reactive  species 
over  the  top  of  the  barrier  or  tunnelling  through  the  barrier.  The  average 
frequency  <u>  reflects  this  situation.  Because  the  Metropolis  sampling 
method  deals  with  discrete  collections  of  particles,  with  specified  forces 
of  interaction,  the  transfer  frequencies  for  over-the-top  of  the  barrier  and 
tunnel  transfers  can  be  determined  in  terms  of  the  actual  interactions  used 
instead  of  using  non-specific,  model  potential  energy  functions  for  the 


barrier . 


Introduction 


Numerical  saddle  point  methods  are  available  to  determine  the  location 

of  the  transition  state  for  chemical  reactions  which  involve  atom 

transfer;1  these  methods  generally  have  been  used  to  determine  the 

configuration  of  the  transition  state  as  a  single  configuration  on  an 

adiabatic  potential  energy  surface  for  a  reaction  in  vacuum.  Although  these 

calculations  are  routinely  carried  out,  they  have  not  been  generally  applied 

to  reactions  which  take  place  in  solution.  Our  purpose  is  to  develop  an 

algorithm  which  combines  saddle  point  optimization  with  Metropolis/Monte 
2  3 

Carlo  sampling  ‘  in  order  to  determine  values  of  rate  constants  for 
reactions  in  solution.  No  new  theory  is  developed  here;  we  work  within  the 
structure  of  the  transition  state  theory.  tfe  assume  that  the  relatively 
simple  Arrhenius  form  of  the  rate  constant, 

k  -  A  exp(-EykT).  (1) 

applies  and  that  it  is  possible  to  estimate  the  activation  energy  E  as  the 
difference  of  the  average  energies  of  the  transition  state,  <E^>,  and  the 
initial  state  <EQ>.  The  pre- exponential  frequency  factor  is 
predominantly  an  average  of  the  frequency  of  passage  through  the  transition 
state . 

An  important  aspect  of  this  work,  we  believe,  is  the  fact  that  it  is 
possible  to  calculate  the  activation  energy  and  the  pre -exponential 
frequency  factor  in  terms  of  the  discrete  configurations  of  particles 
generated  in  the  Metropolis  sample  together  with  the  actual  forces  which  one 
chooses  to  represent  the  interparticle  Interactions.  It  is  not  necessary, 
for  example,  to  use  non-specific  (and  often  one -dimensional)  models  of  the 
barrier  in  the  transition  state.  As  a  result,  it  is  possible  to  carry  out 


calculations  which  are  more  sensitive  to  the  structure  and  nature  of  the 
interactions  in  the  reactive  system  plus  solvent. 

In  the  following  sections,  we  discuss  the  evolution  of  a  reaction  with 
the  use  of  the  Monte  Carlo  simulation.  An  important  issue  is  the  weighting 
of  the  distribution  in  a  way  which  samples  the  transition  state.  Given 
adequate  sampling,  it  is  possible  to  examine  each  configuration  to  determine 
the  how  the  reactive  system  passes  through  the  transition  state:  either  by 
passing  over  the  barrier  or  by  tunnelling  through  it.  The  methods  are 
illustrated  with  a  simple  example,  namely,  the  inversion  of  a  bent  molecule 
of  the  form  ABA  in  a  two-dimensional  solvent  of  Lennard- Jones  disks. 


Metropolis  Sampling  and  the  Transition  State 


Metropolis  sampling  is  an  energy  minimization  technique.  This  is  an 
important  consideration  when  one  attempts  to  sample  configurations  of  the 
transition  state  ensemble.  The  Metropolis  sampling  algorithm  is  well  known 
and  documented.  Averages  of  various  quantities  of  interest  can  be  computed 
directly  as  the  ensemble  of  systems  evolves: 

1  " 

<q>  -  5  l  \  (2) 

i-i 

for  N  samples  where  q^  is  the  value  of  the  variable  q  in  the  1th 
configuration  of  the  Metropolis  sample. 

There  are  two  ways  to  apply  Metropolis  sampling  to  the  investigation  of 
the  transition  state.  In  the  first  method,  a  saddle  point  optimization  to 
the  transition  state  is  calculated  from  the  Initial  state  configuration. 

This  is  clearly  the  most  rigorously  correct  approach.  However,  it  is  the 
most  costly  path;  saddle  point  optimizations  are  not  generally  fast  or 
efficient.  On  the  other  hand,  for  the  second  method,  the  distribution  can 
be  biased  so  that  configurations  of  the  transition  state  are  sampled  in  the 
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same  manner  as  one  samples  configurations  of  the  initial  state.  This 
approach  is  computationally  much  more  efficient.  The  latter  approach  is  the 
one  we  take. 

Ve  consider  reactions  which  can  be  followed  by  following  the  migration 
of  a  definite  atomic  or  molecular  species.  The  inversion  of  a  bent  molecule 
ABA  in  a  two-dimensional  solvent  is  one  example: 


(a) 


Another  example  is  the  transfer  of  an  atom  in  a  bimolecular  reaction  of 
the  type 

X-Y  +  Z  — >  X  +  Y-Z  (b) 

In  each  ca.e,  the  location  of  a  particular  atom  with  respect  to  the  others 
allows  one  to  follow  the  reaction;  the  atom  B  in  the  ABA  inversion  and  the 
atom  Y  in  the  XY  +  Z  to  X  +  YZ  transfer  will  be  referred  to  as  "transfer 
species"  in  the  remainder  of  the  paper. 

Given  any  individual  configuration  of  the  transition  state  as  generated 
in  an  appropriate  Metropolis  cycle,  the  transfer  species  can  execute  only 
one  of  two  types  of  motion  in  the  transition  state:4  locally  stable  or 
locally  unstable  motion.  The  stability  of  the  motion  is  measured  along  to 
the  transfer  axis  and  is  determined  with  respect  to  the  interaction  of  the 
transfer  species  with  all  the  remaining  particles  in  the  instantaneously 
rigid  surroundings  of  the  configuration. 

In  the  first  case,  the  transfer  species  occupies  a  local  energy  well  in 
the  Instantaneous  configuration  of  all  the  other  atoms  in  the  transition 
state.  Ve  assume  that  a  configuration  in  which  the  transfer  species  is 
carried  to  and  on  through  the  transition  state  in  a  local  energy  well 
contributes  to  the  totally  adiabatic  limit  of  the  reaction;  the  transfer  is 
both  electronically  and,  with  respect  to  the  solvent,  vibratlonally 


adiabatic.  On  the  other  hand,  if  the  transfer  species  is  locally  unstable 
in  the  transition  state,  it  must  actually  reach  some  state  near  to  the 
transition  state  and  tunnel  through  the  remainder  of  the  barrier. 


In  view  of  the  fact  that  there  is  only  one  of  two  possible  modes  of 
transfer  through  the  transition  state  for  each  configuration,  the  jump  rate 
v  can  be  expressed  as  the  simple  Arrhenius  formula 

v  -  <w>exp[-0<<E*>-<Eo>)]  (3) 

in  which  the  energy  difference  is  the  activation  energy  and  <u>  is  the 
average  value  of  the  frequency  factor  for  the  transfer,  -  1/kT: 


<w> 
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w  for  an  "over  the  barrier"  transfer 

« 

w  for  a  tunnel  transfer, 
d 


(4) 


(5) 


In  the  next  section,  we  outline  the  test  for  local  mechanical  stability 
of  the  transfer  species.  Following  that  discussion,  we  illustrate  the 
technique  by  applying  it  to  the  case  of  the  inversion  of  SOz  for  which  there 
is  a  good  analytical  potential  energy  function.  In  the  remaining  sections 
of  the  paper  we  consider  the  calculation  of  the  transfer  frequency  for 
passage  over  the  barrier  and  for  tunnelling  through  the  residual  barrier. 
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Analysis  of  Mechanical  Stability 


The  transition  state  is  characterized  by  global  mechanical  instability. 
Nevertheless,  it  is  possible  for  the  migratory  atomic  or  molecular  species 
to  occupy  a  state  in  the  transition  state  which  is  locally  mechanically 
stable.4  It  is  necessary  to  be  able  to  test  any  given  configuration  of  the 
ensemble  of  transition  states  to  see  if  this  condition  is  satisfied. 

It  is  possible  to  consider  the  stability  of  the  motions  of  the  transfer 
species  within  a  given  field  of  the  surrounding  atoms  of  the  reactants  in 
their  configuration  in  the  transition  state.  The  determiniation  is  made 
using  the  second  order  coefficient  of  the  Taylor  series  of  the  potential 
energy  function  of  the  system  expressed  in  terms  of  symmetry  adaptable 
functions,  namely,  the  spherical  harmonic  functions.5 

The  interactions  between  atoms  within  a  molecule  can  be  expressed 
adequately  in  terms  of  a  collection  of  two-center,  three-center  and  higher 
order  functions.  The  two-center  contributions  to  the  potential  energy 
function  can  assume  familiar  forms  such  as  the  Morse,  Lennard- Jones  or 
Rydberg  functions.  The  three-center  and  higher  order  terms  can  be  expressed 
as  polynomial  functions  in  terms  of  the  bond  variables ,  p ^  -  R°  and 

a  range  function  that  ensures  the  appropriate  limiting  behavior  on 
dissociation  (cf.  Murrell,  ref.  6). 

A  general  potential  energy  function  can  have  an  angle - dependence : 5 


F(r)  -  f(r)Y^<r) 


(6) 


4 


in  which  the  angular  part  is  contained  in  the  spherical  harmonic  function 

A 

Y^(r)  and  f(r)  is  the  radial  part.  Complicated  angular  dependencies 


clearly  can  be  constructed  as  sums  of  products  of  spherical  harmonic  and 
radial  components.  A  scalar  function  has  the  simple  representation 
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V(r)  -  ^5irY  (r)v(r) 
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The  expansion  of  Chis  particular  (scalar)  function  as  a  Taylor  series  has 


the  form 


®>  n 

k3/2r  V  /n\  ““*/  - 


V<R+Vri>  -  (to)"*  I  I  (•)  ^ 


n“0  i ■ 0 


X  [  (2*^1)  (2£  +1)  <2£*-l) ]'1/a(£  £  ■  m  I  fa) 

{1, a)  '  4 
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for  which  A  -  is  given  by 


A  .  - 


n! (2£+l) 


n l  (n-£) ! ! (n+f-1) ! ! 


for  n  £  l  and  n  -  l  - 


for  n  <  l  and  n  -  l  -  odd 


and  I^R)  is  given  by 


‘m<R>  -  ‘-^‘(ns)  s  [ai)n  *R  V<R> 


(it  ■  a  |  fa)  is  the  Clebsch-Gordan  coefficient  and  (5)  is  the  binomial 
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coefficient. 


At  any  given  point  of  expansion,  it  is  possible  to  transform  to  a 


Cartesian  representation.  For  the  first  and  second  order  terms,  in 


particular,  one  finds 


VVi'  -  (r  -  r  )*R  I  (R) 
1  2  11 


where  R  is  the  unit  vector,  and 


V(2)  -  (r  -  r  )TK(r  -  r  ). 
1  2  1  2 


In  a  Cartesian  representation,  the  elements  of  the  force  constant  matrix  K 
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for  the  Cartesian  diagonal  terms  and 
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k  -  X  X  I  (R)/R 

ij  1  J  22  ' 


(14) 


for  the  Cartesian  non-diagonal  contributions.  In  these  expressions,  R  is 


the  distance  between  the  points  I  and  J.  The  lower  case  indices  cover  the 
three  Cartesian  coordinates  x,  y,  and  z. 

If  I  -  0  is  the  label  of  the  transfer  species  within  the  field  of  the 
atoms  in  the  transition  state,  then  stability  of  that  species  is 
determined  by  the  examination  of  the  eigenvalues  of  the  3X3  matrix  of  local 
force  constants.  The  elements  of  the  stability  matrix  are  found  with 


Let  the  eigenvalues  of  the  k-matrix  be  where  i-1,2,3.  Stable  motion  of 
the  transfer  species  in  the  transition  state  is  guaranteed  only  if 

#c  >  0,  all  i  (16) 

In  this  case,  then,  the  transfer  species  occupies  an  energy  well  in  the 
transition  state  configuration.  The  entire  transition  state,  however,  is 
still  globally  unstable. 

On  the  other  hand,  if  <  0  for  any  single  i  (-1,2,3),  the  transfer 
species  cannot  occupy  an  energy  well  in  the  transition  state.  It  is 
necessary  for  the  particle  to  tunnel  through  the  remaining  barrier  to  the 
final  state.  This  problem  is  considered  shortly. 

Example:  Stability  of  the  Linear  Transition  State  of  Sulfur  Dioxide 

In  order  to  illustrate  the  test  for  mechanical  stability,  we  apply  it 
to  a  simple  inversion  in  sulfur  dioxide  for  which  a  good  potential  energy 
function  exists,  namely  the  Murrell  ec  al.  potential.8,7 

Consider  the  configurations  illustrate  in  Figure  1.  The  initial  state 


of  S(>2  is  benC.  In  vacuum,  the  transition  state  for  the  inversion  is 
linear.  By  symmetry,  the  motion  of  sulfur  along  the  0-0  line  will  be 
stable.  The  motion  of  sulfur  along  the  transfer  axis  perpendicular  to  the 
0-0  line  may  or  may  not  be  stable.  As  a  result  of  the  symmetry  of  the 
system,  the  axial  force  constant  can  be  expressed  asS* 


k  -  -  i  I  [I  (R  )  -  2P  (cos  8  )I  (R  )] 
2  3  u  1  20  I  2  HI  22  I  J 


where  P2(cos  8)  is  the  second  order  Legendre  polynomial  and  tf  is  the  angle 
between  the  vector  Riq  and  the  transfer  (z)  axis.  The  quantities  I  and 
I  are  simply8 

T  d2V  .  2  dV 

VR)  '  —2  *  R  dR  (18) 

OK 


i  (R)  -  O  .  I  dv 

a<  W  dR2  RdR  ' 


The  Murrell  potential  is  the  sum  of  contributions.  The  first  terms  are 
pairwise  Rydberg  terms 


V  (R)  -  -  D[l  +  E  a  p1|exp(-a  a) , 
1  i-i  J 


where  D  is  the  dissociation  energy,  p  -  R  -  Rq  and  the  a^  are  coefficients 
in  the  functional  form  of  the  potential.  The  remaining  contributions  come 
from  many-body  terms: 

VR)  -  Vo  TT  [ 1  -  tanh<7  S  /2)]  [l  +  l  c  p  +  £  cop  +  .  .  .1  (21) 

i-i  '■i  iij  J  3  J 

R  is  the  set  of  atomic  coordinates .  The  quantity  is 

S  -  T  b  p  (22) 

l  "  ijKj 

and  the  transformation  matrix  b  is  determined  by  symmetry.  The  remaining 
quantities,  7^  and  c^,  c  ,  etc.  are  determined  by  fitting  the  potential  to 
various  data- -spectroscopic  and  ab  Initio  calculations.8,7 

The  values  of  the  various  coefficients  were  taken  from  Carter  et  al . 1 
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The  R _  distance  is  simply  twice  the  R„  distance  in  the  transition  state. 

00  so 

The  energy  of  the  system  in  its  linear  transition  state  is  plotted  in  Figure 
2  as  a  function  of  the  R^  distance.  The  equilibrium  SO  bond  length  is 
1.431  A.  There  are  two  energy  minima  along  the  0S0  linear  expansion  in  the 
transition  state:  the  first  is  approximately  1.35  A.  Figure  3  shows  the 
stability  of  the  motion  of  the  sulfur  atom  with  reference  to  a  displacement 
along  the  transfer  axis  in  the  transition  state.  It  is  clear  from  the 
figure  that  at  the  lowest  energy  of  the  system  in  the  transition  state, 
approximately  1.35  A,  the  sulfur  atom  cannot  occupy  a  position  of  stable 
mechanical  equilibrium.  This  is  true  for  both  the  pair  potential  and  the 
many-body  contribution.  As  a  result,  one  anticipates  that  any  inversion  of 
the  S02  molecule  would  involve  the  tunnelling  of  sulfur  through  the  0-0 
line. 

An  Interesting  feature  of  Figure  3  is  the  relatively  small  range  of 
values  of  R^  which  admit  over- the -barrier  transfers  of  the  sulfur  atom. 

The  stability  of  the  Rydberg  part  of  the  whole  potential  energy  function  is 
standard  for  pair-potentials;  the  behavior  of  the  Rydberg  potential  is 
similar  to  that  of  the  Morse  or  Lennard- Jones  potentials.  The  many-body 
part  of  the  whole  potential  function  does  not  admit  any  stable  configuration 
for  sulfur  in  the  linear  transition  state.  This  can  be  seen  to  be  the 
result  of  a  strong  restoring  force  for  the  bending  mode.  The  combination  of 
pair-  and  many-body  (angle  bend)  forces,  however,  does  allow  a  small  set  of 
configurations  for  which  it  is  possible  for  sulfur  dioxide  to  invert  by 
carrying  the  sulfur  atom  through  the  transition  region  in  a  potential  energy 
well.  Such  a  state  involves  the  excitation  of  the  SO  stretch  displacements 
to  a  considerable  extent.  As  can  be  seen  in  the  Figures,  there  is  a  second 
minimum  configuration  of  the  transition  state  at  about  1.85  A.  It  is 
unlikely  that  at  ambient  temperatures,  for  example,  there  would  be  any 
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significant  population  of  states  characterized  by  the  lengthened  0-S  bonds. 

Therefore,  on  the  basis  of  this  examination,  we  are  inclined  to  believe  that 

sulfur  would  have  to  tunnel  through  the  0-0  barrier  in  most  cases. 

Similar  stability  conditions  can  be  seen  to  arise  out  of  more 

approximate  representations  of  the  angle  bending  motions.  Smith  and 

Overend,  for  example,  examined  several  potential  energy  functions  as 

representations  of  the  bending  potential  in  bent  triatomic  molecules .  Given 

the  fact  that  the  potential  for  the  inversion  of  S02  must  be  a  symmetric 

double  well  function.  Smith  and  Overend  examined  the  following  form  which  is 

91) 

suggested  by  the  Swalen-Ibers  treatment  of  ammonia: 

V(*-w)  -  k (8-n)2  +  b  exp[ -c(tf-w)2]  (23) 

and  the  angle  8  is  measured  from  the  linear  transition  state  for  the 
inversion.  The  coefficients  k  and  b  depend  on  the  SO  bond  lengths: 
k  -  K  exp[A(Ri+R2)] 

b  -  B  exp[/J(Ri+R2)]  (24) 

and  A  <  0,  etc.  Using  eqs  (23)  and  (24)  in  the  equations  to  determine  the 

stability  in  the  linear  triatomic,  one  finds 

k  -  .xp(2^Rao)  (25) 

so 

[note,  in  the  transition  state,  8  -  w] .  As  0  <  0,  it  is  clear  that  k  <  0  in 
accord  with  the  results  found  above  for  the  Murrell  potential6  (p  in  this 
expression  is  not  1/kT) .  As  a  consequence,  accurately  determined 
experimental  force  constants  for  the  bending  motions  in  a  molecule  can  aid 
in  the  determination  of  the  stability  of  the  transition  state  of  a  reaction. 

Frequency  Factor  for  Completely  Adiabatic  Transfers 

Ve  consider  the  case  for  which  the  transfer  species  occupies  a  stable 
energy  well  in  a  single  transition  state  configuration  of  the  Metropolis 
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sample.  The  problem  is  Co  determine  Che  frequency  factor  which  is 
associated  with  the  migrating  particle  in  this  case. 

The  test  for  stability  of  the  migratory  transfer  species  requires  the 
determination  of  the  eigenvalues  of  the  iorce  constant  matrix  for  the  motion 
of  the  transfer  species  in  the  fixed  field  of  the  surrounding  atoms  of 
solute  and  solute.  If  all  of  the  eigenvalues  are  positive  definite,  all 
motions  of  the  transfer  species  are  locally  mechanically  stable.  The 
oscillatory  motion  of  the  transfer  species  along  the  axis  in  the  direction 
of  the  final  state  determines  the  frequency  which  is  associated  with  the 
wholly  adiabatic  transfer.  Let  the  z-axis  be  aligned  with  the  tangent  to 
the  transfer >axis  in  the  transition  state.  The  force  constant  for  harmonic 
motion  of  the  migrating  atom  along  this  axis  is 

K  -  1  l  <I20<R  >  +  2P2(cos  «ri>I22<V>  (26) 

in  which  the  summation  is  taken  over  all  atoms  in  the  surroundings.  The 
frequency  for  the  transfer  is  simply 

-  2VR7S,  (27) 

where  the  factor  of  2  arises  via  the  same  argument  as  used  in  the  method  of 
beats,10  and  m  is  the  mass  of  the  transfer  species.  This  frequency  is  used 
for  a  particular  configuration  in  eq  (5). 


Frequency  Factor  for  Tunnel  Transfers 


The  accuracy  of  the  frequency  factor  for  a  nuclear  tunnel  transfer 
depends  on  the  accuracy  with  which  one  can  describe  the  system  in  the  first 
place.'  Extensive  use  has  been  made  of  one  dimensional  model  potential 
energy  functions  as  representations  of  the  barrier.  We  have  been  able  to 
show11  recently  that  it  is  possible  to  get  good  agreement  with  experiment 
for  the  inversion  tunnelling  in  the  ammonia  by  using  a  representation  which 
employs  discrete  pair -potentials  for  the  NH  interactions.  Thus,  in 
principle,  if  good  model  potential  energy  functions  are  known  for  the  atomic 
pair  interactions  and  the  angle  bending  modes,  it  is  possible  to  carry  out 
accurate  calculations  to  estimate  the  tunnel  frequency. 

We  propose  the  following  computational  algorithm.  Given  a  case  for 
which  tunnelling  is  required,  locate  the  nearest  positions  of  local 
stability  in  the  regions  of  the  initial  and  final  states  for  the  transfer 
species  in  the  frozen  configuration  of  the  environment.  [It  is  necessary  to 
note  that  the  initial  and  final  states  for  the  tunnel  transfer  are  almost 
never  the  same  as  the  thermodynamic  initial  and  final  states  of  the 
reaction. ]  The  rationalization  we  use  is  that  a  fluctuation  of  the  system 
will  bring  the  transfer  species  close  to  the  transition  state.  These 
positions  of  local  mechanical  stability  can  be  found,  for  example,  with  the 
use  of  steepest  descents  routines. 

Having  determined  stationary  initial  and  final  states  for  the 
calculation,  it  is  necessary  to  carry  out  some  form  of  separation  of  the 
local  modes  from  those  of  the  remainder  of  the  system.  It  is  also  necessary 
to  identify  the  perturbations  which  account  for  the  tunnel  transfer. 

We  assume  that  Born-Oppenheimer  electronic  states  have  been  separated 
from  molecular  vibrations;  the  reactive  subsystem  migrates  over  an  adiabatic 


potential  energy  surface.  Ve  identify  and  separate  local  vibrational  modes4 

of  the  reactive  solute  from  the  remaining  solvent  with  a  second  use  of  the 

12 

Bom-Oppenheimer-Holstein  adiabatic  separation  applied  to  the  intra-  and 
interaolecular  vibrations.  The  procedure  we  follow  for  atomic  motions  is 
analogous  to  the  method  which  is  used  for  the  electron  transfer 
problem. 13,14  The  method  works  for  the  atom  transfer  for  the  reason  that 
atom- tunnel ling  only  occurs  if  the  vibrational  overlap  factors  are 
nonvani shingly  small  and  interaction  with  solvent  is  relatively  weak.  This 
is  so,  of  course,  only  for  reasonably  small  separations  between  the  initial 
and  final  states  for  the  migration. 

The  reactive  subsystem  can  be  defined  to  consist  of  the  underlying 
reactive  solute  and  nearest  neighbor  solvent  when  that  solvent  Interacts 
strongly  with  the  solute.  The  point  we  make  here  is  that  molecular 
vibrations  can  be  handled  (certainly  in  the  harmonic  limit)  reasonably 
accurately,  if  the  number  of  atoms  and  molecules  is  not  too  large.  Thus,  it 
is  possible,  in  principle,  to  define  the  reactive  subsystem  to  consist  of 
solute  and  sufficient  solvent  such  that  were  any  additional  solvent 
considered,  its  effect  would  merely  mimic  the  bulk  properties  of  the 
solvent. 

A  general  development  of  the  Bom-Oppenheimer-Holstein  separation  of 
local  vibrational  modes  for  a  reactive  system  has  been  presented 
elsewhere;4*  we  summarize  that  development.  A  plane  is  located  at  the 
saddle  point  for  the  reaction.  Stationary  initial  and  final  states  are 
located  on  either  side  of  this  plane.  Functions  associated  with  these 
states  are  basis  functions  for  the  description  of  the  tunnel  transfer.  It 
is  possible,  as  is  shown  in  the  Appendix,  to  include  basis  functions  of  the 
transition  state  itself  in  order  to  account  for  tunnel  transfers  which  are 
oblique;  that  is,  transfers  which  do  not  occur  on  a  straight  line  from  the 


initial  to  th«  final  state  through  the  transition  state. 


The  Hamiltonian  operator  for  the  complete  system  of  reactant  and 
surroundings  is  written  simply  as 

H  -  H  +  H  +  V"  (28) 

f  n  r 

in  which  is  the  vibrational  Hamiltonian  operator  for  the  reactant 
subsystem  and  H  is  the  operator  for  the  remainder  of  the  system.  The 

•n 

specific  interactions  between  individual  atoms  of  the  reactive  subsystem  and 
atoms  and  molecules  of  the  surrounding  environment  are  contained  in  the  term 
V*11;  these  individual  contributions  can  be  enumerated. 

t 

As  an  example  of  the  Born-Oppenheimer-Holstein  separation,  we  consider 
the  specific  separation  for  the  two-state  system  which  applies  when  the  the 
transfer  takes  place  only  between  ground  initial  and  final  vibrational 
states.  The  specific  formulae  that  result  should  be  generally  applicable 
because  of  the  great  probability  of  a  tunnel  transfer  between  the  ground 
initial  and  final  states  which  have  been  adiabatically  formed.  Locally,  one 
finds 


Hr*I(F)  “  eI(F)*I(F)  (29) 

In  terms  of  a  basis  set  which  is  orthonormal  with  respect  to  the  initial  and 

final  states,  the  overlap  matrix  is 


S  - 


(30) 


The  inverse  is 


S 


-1 


1 


(31) 


The  state  function  is  constructed  as 

*  -  fjXj  +  (32) 

The  expansion  coefficients  f  are  functions  of  the  set  of  environmental 
coordinates  of  the  solvent  and  molecular  framework  of  the  reactive  solute. 
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In  Ch«  usual  unntr,  writ* 

Ht  -  E* 

and 

<Xt|H|*>  -  E<Xj | *> 

Upon  expansion. 


(33) 


(34) 


E(f  +  S  f  )  -  <x  |H  +  V*“|X  >f  +  <x  |H  +  V*"|x  >f 

I  OJ  j'  AI*  r  ■*!  I  *I>  r  |AJ  J 


-<*.!»„ U,>£:  <35> 

V*“  contains  all  eh*  interactions  which  depend  upon  coordinates  of  the 
reactive  subsystem  and  the  environment.  Rewrite  eq.  (35)  as 

IS  (E  -  <  -  H  )f  -  XlV  (1  -  S  )  +  V,n)f  (36) 

J  ij  i  an  1  **  1J  ij  ij  i  v  ' 

This  equation  can  be  written  more  transparently  as  the  matrix  equation 

SHF  -  VF  (37) 

where  H  is  the  matrix 


H 


fE-H  - c  0  \ 

0  E-H  -«J 

•n  hr 


(38) 


In  this  equation,  is  the  Hamiltonian  operator  for  the  solvent  species 
without  the  specific  interaction  between  solute  and  solvent.  F  is  the 
vector 

ft 


F  - 


f 

s  b. 


(39) 


and  V  is  the  matrix 


V  - 


V*“  V  +V’ 


V  +v 

ba  ba 


ab  ab 

.•a 


bb 


(40) 


In  this  expression,  V  without  the  superscrip  index  "en"  indicates  a  matrix 


elesMnt  which  involves  the  interaction  between  the  initial  and  final  state 
configurations  and  the  migratory  species.  These  interactions  are  contained 


a 


a 
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within  the  operator  H^.  We  have  Ignored  matrix  elements  which  involve  the 
momentum  operators.  They  can  easily  be  Included  in  this  general  term.  In 
analogy  to  the  electron  transfer  case,  the  matrix  elements  V  account  in 
large  part  for  the  perturbations  which  drive  the  tunnel  transfer. 


At  this  point,  one  can  write 
HF  -  S^VF 


(41) 


Upon  expansion  and  rearrangement  of  this  equation,  we  find  the  equation  of 
motion  for  the  coefficient  f  ,  for  example, 


fE  -  H  -  *  -  - (V*B  -  S(Vv  +  V*n))lf 

«n  a  ^  j2  ba  ba  j  a 


- i — f 

1  -  S2^ 


V  +  v"  -  svHf 

ab  ab  bbj  b 


(42) 


With  the  use  of  this  equation,  the  matrix  element  Lip  which  accounts  for  the 
tunnel  transfer  can  be  identified: 


-  -—M 

•b  1  -  S2  l 


V  +  V  -  sv 

ab  ab  bb 


(43) 


with  a  similar  term  for  L  for  the  reverse  transfer. 

ba 

For  an  individual  configuration  of  the  ensemble,  a  transfer  across  the 

o 

barrier  takes  place  in  the  0  K  limit.  The  tunnel  frequency  is  calculated 
with  the  use  of  the  "Golden  Rule"  expression: 


2  w 


w  -  ^  Z  |L  |«(E  -E  ) . 

if  ft  u  1  if'  i  f 


(44) 


Assuming  harmonic  environmental  modes  and  the  Condon  approximation  for  the 
matrix  element  L^,  it  is  possible  to  complete  the  summation  over  final 


states.  The  result  is 
2* 


/  2  2  v  d 

(a  q„/2) 

ftwn! 


where  u  is  the  frequency  of  the  accepting/donating  mode  and 

a  -  /mw/ft 


(45) 


(45) 


i 

and  qg  is  Che  (dimensionless)  equilibrium  coordinate  of  Che  normal 
accepting/donating  mode.  The  number  n  is  fixed  by  energy  conservation: 

| Eg |  -  nAw  (46) 

and  Eq  is  the  energy  difference  between  the  initial  and  final  states. 

The  energy  Aw  in  the  denominator  of  eq  (45)  comes  from  the 
approximation  of  the  density  of  states  of  the  manifold  as  1 /Aw  (Siebrand, 
ref  15).  This  is  an  approximation  which  has  proved  to  be  inadequate  for  the 
radiationless  transition  in  benzene.  It  is  clear  that  some  better 
approximations  will  have  to  be  used  to  obtain  accurate  estimates  of  the 
tunnel  frequency.  For  the  moment,  the  expression  (45)  serves  to  identify 
the  route  which  we  suggest  to  follow  in  order  to  determine  the  transfer 
frequency.  There  is  indeed  another,  even  more  approximate,  route  to  follow: 
one  can  use  the  simple  method  of  beats  for  which  the  transfer  frequency  is 
given  by10 


w  - 

if 


2ILJ 


(47) 


This  expression  will  be  reasonably  accurate  whenever  the  energy  difference 
between  the  initial  and  final  state  energy  wells  is  vanishingly  small. 

This  expression  is  the  one  which  we  use  in  the  model  simulation  of  the  next 
section. 


Example:  Rate  of  Inversion  of  a  Bent  ABA  Molecule 

The  reaction  we  simulate  is  the  hypothetical  inversion  of  a  bent  ABA 
molecule,  see  reaction  (a).  The  mass  of  the  B-species  is  taken  arbitrarily 
to  be  12  au,  viz.,  carbon.  The  purpose  of  this  choice  is  to  try  to 
determine  whether  an  atom  of  the  mass  of  carbon  tunnels  to  any  extent. 

The  following  interactions  are  used  between  the  atoms  of  the  solute  and 


between  the  solute  and  solvent.  First,  the  Lennard- Jones  interaction  is 


used  Co  describe  Che  bonds  beCween  A  and  B.  The  inceraccion  between  Che  two 
atoms  A,  however,  is  modelled  as  harmonic.  The  reason  for  doing  so  lies 
with  Che  fact  Chat  a  Lennard- Jones  potential  is  Coo  "soft"  in  the 
dissociative  region  Co  account  for  Che  effect  of  other  atoms  of  solute  and 
solvent  which  in  a  real  system  would  present  an  effective  repulsion  to  the 
expansion  of  the  A-A  distance.  Thus,  with  the  use  of  the  harmonic 
interaction,  as  opposed  to  the  Lennard-Jones  which  would  be  more  appropriate 
to  the  species  in  the  gas  phase,  it  is  possible  to  model  the  effect  of  an 
atom  transfer  within  a  rigidly  bound  polyatomic  molecule  with  the  use  of  a 
minimal  system  as  we  have  done  here.  Note,  we  have  not  used  angle  bend  or 
higher  order  interactions  in  this  calculation,  as  correctly  we  should  do. 

The  reason  for  omitting  them  is  primarily  due  to  the  fact  that  we  have  not 
yet  worked  out  a  good  way  to  approximate  the  interaction  in  a  form  which  can 
be  incorporated  into  the  tunnel  calculation.  We  are  currently  seeking  ways 
to  approximate  these  important  interactions. 

The  solvent  is  modelled  as  a  collection  of  uniform  spheres  of  radius 
2.8  A.  A  Lennard- Jones  potential  is  also  used  to  describe  the  interaction 
between  the  inidividual  atoms  of  the  solute  and  the  solvent.  Finally,  the 
interactions  between  the  molecules  of  solvent  are  also  modelled  as 
Lennard- Jones.  The  parameters  are  summarized  in  Table  1. 

We  used  Lennard-Jones  potentials  throughout  the  calculation  in  the 
interest  of  computational  speed.  The  form  of  the  potential  used  was 

Vu  -  « (a/r)*£(a/r)8  -  2]  (48) 

in  which  a  is  the  distance  for  which  -  -e.  More  elaborate  potential 
energy  functions  usually  require  calls  to  exponential,  error  and  similar 
functions.  This  fact  can  greatly  increase  the  amount  of  time  needed  to 
carry  out  a  large  number  of  calculations . 

In  carrying  out  the  calculation,  we  fixed  the  central  B  atom  in  the  ABA 
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molecule  at  Che  center  of  Che  unit  cell  In  Che  Metropolis  system.  The 
A-aComs  and  Che  solvent  were  free  to  move.  The  entire  system  was 
two-dimensional;  the  solvent  was  modelled  as  a  collection  of  Lennard- Jones 
disks.  In  order  to  eliminate  edge  effects,  the  minimum  image  technique  was 
used.  The  unit  cell  was  surrounded  by  eight  image  cells.  Each  particle 
within  Che  primary  cell  interacted  either  with  another  species  in  the  cell 
or  its  nearest  image  if  the  image  was  closer.  Sampling  was  carried  out  on  a 
canonical  ensemble. 

The  calculations  required  the  use  of  the  Cray  XMP/24  computer  in  order 
to  obtain  results  in  reasonable  time.  The  Metropolis  simulation  proceeded 
in  two  steps.  First,  the  unconstrained  optimization  of  the  system  was 
carried  out  in  order  to  determine  the  average  energy  of  the  initial  state . 

No  other  averages  were  determined.  These  optimizations  were  carried  out  for 
a  sequence  of  temperatures  beginning  with  a  cycle  at  4  K  to  300  K.  The 
simulation  of  the  transition  state  involved  calculations  for  the  same 
sequence  of  temperatures  as  used  in  the  determination  of  the  initial  state. 
The  algorithm  to  determine  the  energy  and  frequency  factor  for  the 
transition  state,  however,  was  considerably  more  complicated  than  that  used 
in  the  determination  of  the  initial  state. 

The  vacuum  transition  state  for  the  ABA  system  is  the  optimized  linear 
configuration  A-B-A.  With  the  B- species  fixed  to  the  center  of  the  unit 
cell,  and  the  A-species  constrained  to  move  along  the  A-B-A  line,  the  system 
was  allowed  to  evolve  configurations  of  the  transition  state  in  the 
following  manner.  Motion  of  the  A-species  along  the  A-B-A  line  was  allowed 
together  with  motion  of  the  solvent  to  optimum  configurations.  We  made  the 
assumption  that  the  solute  in  its  transition  state  could  adiabatically 
adjust  to  the  solvent.  Thus,  the  Metropolis  steps  generated  configurations 
of  solute  and  solvent  which  were  near  to,  but  not  exactly,  the  transition 


state.  In  order  to  generate  the  final  configuration  which  represented  the 
transition  state,  we  carried  out  a  Newton -Raphs on  saddle  point  optimization 
of  the  ASA  solute  in  the  frozen  solvent.  Ue  reasoned  that  the  solvent  would 
equilibrate  to  a  configuration  close  to  the  transition  state.  Therefore, 
the  equilibration  of  the  solute  in  the  rigid  solvent  would  not  return  the 
solute  to  a  ground  or  initial  state  configuration.  In  carrying  out  this 
kind  of  optimization  to  the  transition  state,  it  is  clear  that  geometrically 
non-linear  transition  states  are  possible. 

For  each  new  configuration  of  solvent  generated  in  the  Metropolis 
sample,  a  Newton- Raphson  optimization  was  carried  out  for  the  B- species 
alone  to  determine  the  nearest  stationary  states  on  the  reactant  and  product 
sides  in  the  frozen  configuration  of  the  transition  state.  Note,  these 
initial  and  final  states  are  not  the  thermodynamic  initial  and  final  states. 
They  are  states  below  the  transition  state  barrier  in  which  the  B- species 
occupies  a  local  energy  well. 

We  found  that  on  occasion,  the  routine  would  find  two  energy  minima  on 
the  same  side  of  the  A-A  line.  In  this  case,  the  configuration  was 
discarded  and  the  last  configuration  kept.  The  test  of  whether  the  initial 
and  final  states  are  separated  by  the  A-A  line  is  simple.  The  directed 
distance  from  the  line18 

Ax  +  By  +  C  -  0  (49) 


to  the  point  (x^y  )  is 


d 


1 


Ax  +  By  +  C 

l  J  l _ 

±(A2  +  B2)1'2 


(50) 


In  the  same  manner,  the  distance  to  a  second  point  (x^y^)  from  the  same 

line  is  d  .  If,  then 
2 

S  -  (Axx  +  Byj  +  C)(Ax2  +  Byz  +  C)  >  0  (51) 


both  the  initial  and  final  state  lie  on  the  same  side  of  the  A-A  line.  This 


test  is  easily  generalized  Co  three  dimensions;  Che  A-A  line  is  simply 
replaced  by  the  dividing  plane  which  contains  the  saddle  point  or  some  other 
suitably  defined  point  of  reference. 

The  mechanical  stability  of  the  migratory  B- species  was  determined  for 
each  Metropolis  element  of  the  ensemble  of  transition  states.  If  the 
species  B  proved  to  be  stable  in  the  presence  of  its  environment  of  the 
A-species  and  the  solvent,  that  element  of  the  transition  state  ensemble  was 
adiabatically  attained.  Ue  then  determined  the  frequency  associated  with 
the  motion  of  B  along  the  transfer  axis  according  to  eqs  (26)  and  (27). 

This  frequency  was  identified  as  and  used  in  eq  (4).  If  the  B-species 
was  found  to  be  locally  unstable  in  the  transition  state,  then  a  tunnel 
frequency  was  calculated.  This  frequency  was  determined  by  the  following 
sequence  of  steps. 

The  determination  of  the  tunnel -transfer  frequency  was  found  by  using  a 
one -dimensional  analysis  along  the  transfer- axis .  Thus,  we  considered  only 
the  direct,  straight-line  transfer  from  the  initial  to  the  final  state.  The 
basis  functions  are  the  simple  ground  state  one -dimensional  harmonic 
oscillator  functions: 

uQ(a,z)  -  (*/a)1Mexp(-l/2  az2)  (52) 

and  a  similar  function  for  the  final  state  (identified  in  the  exponent  by 
the  coefficient  fi) . 

Let  z  be  the  transfer-axis.  The  potential  energy  operator  can  be 
10 

witten  as 

CO 

V(z)  -  l  znC  (R)  (53) 

n 

n«0 


where 


In  these  expressions,  R  represents  the  set  of  distances  from  the 
transfer-species  to  the  collection  of  N  solvent  molecules  and  remaining 
atoms  of  the  solute.  The  angle  between  the  z-axis  and  the  location  of  the 


species  i  at  a  distance  from  the  coordinate  origin  is  6. 

The  state  function 

z)  -  f  u(a,z)  +  fu(0,z)  (55) 

•  0  b  0 

yields  the  usual  matrix  elements  which  are  needed  to  evaluate  the  frequency 
factor.  The  overlap  S  is 

S  - 

and  R  is  the  distance  between  the  initial  and  final  locations  of  the 
B-species.  We  specifically  need  to  evaluate  the  matrix  element  V  This 
is 


-  <uo(a,z)|V(z)|uo(0,z)> 

(■  1/2  <*>  « 

-  Vil£E\  exp(-Q2)  l  C  (R  )  l  (S)(-d)n"”[l-(-l)“] 

n-0  °  — n 


m-  0 


i  (m-1) ! !  (.a+fi) 


-m/2 


where 


2  ~2 


Q  - 


and 


d  - 


2(a  +  0) 

al  -  01 

•  b 

a  +  0 


R 


(57) 


(58) 


(59) 


The  quantities  l  and  l  are  the  distances  to  the  initial  and  final 

a  b 

locations  at  a  and  b  as  measured  from  the  origin  on  the  z-axis.  Finally, 
the  quantity  C  (R  )  is  evaluated  with  the  set  of  distances  and  angles 
measured  with  respect  to  the  origin  of  coordinates  on  the  B* species  on  the 
z-axis. 
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„  I 

Ve  made  Che  assumption  Chat  Che  environmental  factors  V  contribute  ; 

only  a  snail  amount  Co  Che  total  value  of  Che  matrix  elemenc  and  may, 
therefore,  be  ignored.  This  assumption  can  always  be  tested  by  direct 

computation  whenever  one  suspects  it  is  not  reasonable.  However,  the  j 

t 

combination  of  small  overlap  and  large  distances  in  C  (R)  would  seem  to 
Justify  the  neglect  of  these  terms  under  most  circumstances. 

The  frequency  of  the  tunnel  jump  was  evaluated  with  eq  (47).  The 
quantities  a  and  were  determined  in  terms  of  the  harmonic  frequencies 
which  could  be  assigned  to  the  transfer- species  in  its  initial  and  final 
locations . 

The  system  consists  initially  of  a  simple  square  lattice  at  4  K. 

Several  runs  were  made  at  4  K  to  "temper*  the  system.  Next,  runs  were  made 
at  50,  100,  200  and  300  K.  This  routine  was  followed  both  for  the 
determination  of  the  initial  state  energy  (unrestricted  Monte  Carlo)  and  for 
the  transition  state  calculations  (biased  in  favor  of  the  linear  transition 
state  for  ABA) . 

The  results  of  the  calculations  are  summarized  in  Table  II.  As  can  be 

seen,  two  values  of  the  (harmonic)  A-A  interaction  were  used,  the  "hard" 

harmonic  Interaction  (viz.,  8  x  10^  dyn  cm  *")  yields  degrees  of  overall 

adiabaticity  of  reaction  which  are  much  smaller  than  those  found  with  the 

4  -1 

smaller  value  of  the  A-A  interaction  (  8  x  10  dyn  cm  ) .  Because  of  the 
weak  A-A  interaction  in  the  first  case,  it  is  clear  that  the  activation 
energy  for  this  system,  when  it  is  positive  at  higher  temperatures,  is  very 
small.  As  one  would  expect,  with  the  stronger  A-A  Interaction,  the 
activation  energy  is  larger,  reflecting  the  work  needed  to  stretch  the  A-A 
distance  in  the  transition  state. 

There  is  a  temperature -dependence  in  the  frequency  factor,  to  a  degree. 

At  low  temperatures,  below  100  K,  the  frequency  factor  increases  with 
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decreasing  temperature.  At  200  X  and  above  it  appears  to  be  essentially 
constant . 

The  degree  of  overall  (electronic  and  vibrational)  adiabaticity  of  the 
reaction  is  easily  determined  in  the  sampling  cycle  by  assigning  a  value  of 
unity  to  any  configuration  which  admits  an  over -the -barrier  transfer  and 
zero  to  the  cases  for  which  tunnelling  through  the  residual  barrier  is 
required.  As  one  expects,  in  the  system  in  which  there  is  a  strong  A-A 
association,  there  will  be  a  greater  amount  of  tunnelling  through  the 
residual  barrier  in  the  transition  state.  This  is  clear  from  the  fact  that 
even  at  300  X,  the  fraction  of  over-the-barrier  transfers  is  only  about 
0.01. 

A  curious  feature  of  these  calculations  is  the  fact  that  especially  for 
the  "weaker"  force  constant  and  at  very  low  temperatures,  the  linear  ABA 
species  seems  to  be  stabilized  by  the  solvent  as  indicated  by  the  apparent 
negative  activation  energies.  It  may  be  the  case  that  at  the  low 
temperatures  in  particular,  a  very  large  number  of  cycles  needs  to  be  used 
in  order  to  reach  an  accurate  value  of  the  initial  state  energy.  It  was 
clear  from  our  results,  however,  that  at  the  higher  temperatures,  the  system 
clearly  stabilizes  for  the  number  of  cycles  used  (from  156,000  to  312,000). 

Vhether  or  not  the  solvent  stabilization  of  the  linear  configuration  at 
low  temperatures  for  the  weaker  interaction  is  real  needs  to  be  examined  in 
more  detail.  We  did  not  determine  average  configurations  of  the  intial 
state.  Ve  did  keep  a  record  of  the  average  cartesian  coordinates  of  the  ABA 
atoms  in  the  transition  state  calculations.  For  the  case  of  the  weak 
A-A  interaction,  we  found  that  the  average  A-A  separation  was  approximately 
2.8  A  which  is  a  very  small  extention  over  the  equilibrium  value  of  2.75  A. 
However,  with  the  stronger  interaction  between  the  A-atoms,  even  this  slight 
extention  is  reflected  in  a  significant  increase  in  energy.  For  the 


3 


a 


3 


4 


A 


$ 


M 

71 


stronger  Interaction,  the  average  A- A  distance  in  the  transition  state  is 
approximately  2.78  A.  For  both  the  "weak"  and  the  "strong"  A-A  harmonic 
interactions,  the  transition  state  remained,  on  average,  linear  in  the 
presence  of  solvent. 

Discussion 

The  objective  of  this  work  is  to  model  the  transition  state  by  a 
combination  of  saddle  point  and  Monte  Carlo  optimizations.  Metropolis 
sampling  generates  individual  elements  of  the  ensemble  of  configurations  of 
the  transition  state.  Saddle  point  optimization  can  then  be  applied  to  each 
configuration  in  order  to  generate  an  ensemble  of  individual  transition 
states. 

The  general  analysis  we  have  outlined  highlights  the  important  features 
of  the  treatment.  Many  issues  need  to  be  explored  and  developed  further. 

An  important  issue,  which  we  have  mentioned,  but  not  included  in  the  example 
calculation,  has  to  do  with  the  role  of  many-body  interactions.  These  are 
arise  predominately  as  angle -bend  interactions  in  vibrational  spectroscopy 
in  the  harmonic  limit.  Bend-bend  and  bend- stretch  interactions  are 
manifestations  of  higher  order  many-body  forces.  The  example  of  SO^,  given 
earlier,  shows  that  these  higher  order  contributions  cannot  be  neglected  if 
one  is  to  get  an  accurate  estimate  of  the  pre- exponential  factor  for  the 
reaction.  It  does  seem  possible,  on  the  basis  of  the  Smith-Overend  work  on 

g 

enharmonic  interactions ,  to  estimate  these  interactions  with  the  use  of 
vibrational  spectroscopy  of  the  initial  and  final  states.  The  systematic 
investigation  of  this  approach  needs  to  be  examined  further,  and  we  plan  to 
do  so. 

There  has  been  other  work  to  apply  Monte  Carlo  methods  to  the  study  of 
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chemical  reactions.  Doll  and  Adams17  in  a  series  of  papers  have 
investigated  the  use  of  Metropolis  sampling  together  with  an  improved 
version  of  Slater's  theory  of  reactivity.  They  have  applied  their 
analysis  to  the  thermal  desorption  problem.  At  first  glance,  our  approach 
is  similar  to  the  Adams-Doll  work.  On  closer  examination,  however,  there 
are  significant  differences.  Doll17*  proposed  examining  the  expression 

L(q)  -  i  <S(x  -  q)x>  (60) 

which  yields  the  average  number  of  zeros  for  the  function  x(t)  -  q  in  an 
interval  of  time  0  £  t  £  r.  The  quantity  q,  in  Slater's  theory,  is  a 
critical  distance  for  reactivity,  a  critical  bond  length  or  height  of  the 
barrier.  Adams  and  Doll17b  d  replace  the  delta -function  by  a  Gaussian  and 
directly  evaluate  the  average  by  the  Metropolis  technique.  The  value  of  q 
is  specified.  The  specification  of  q  is  a  restriction  which  is  not 
satisfactory  for  reactions  taking  place  in  solution.  As  noted  in  the  early 
work  on  the  transition  state  theory,19,20  the  quantity  q  in  essence  is  the 
location  of  the  hyperplane  which  passes  through  the  transition  state  in  the 
reaction  hypersurface.  This  identification  alone,  however,  is  not 
sufficient  to  derive  the  final  form  of  the  transition  state  rate  constant. 
Optimization  needs  to  be  carried  out  on  the  hyperplane  in  order  to  reach  the 
true  saddle  point  which  is  the  transition  state.  It  is  this  prescription 
which  we  have  followed  here.  As  we  have  noted,  it  is  possible  to  assume 
that  a  transition  state  which  is  evaluated  for  the  system  in  vacuum  (in  the 
absence  of  solvent)  can  be  used  as  a  template  on  which  to  bias  the 
distribution  to  the  transition  state  for  the  reaction  in  solution. 

Subsequent  relaxation  of  the  system  to  the  true  transition  state  in  solutin 
is  needed.  As  a  consequence,  we  do  not  see  how  Doll's  approach  can  be 
directly  extended  to  apply  to  the  systems  we  consider. 

In  our  work,  we  have  not  directly  considered  librational  motions  or 
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full  roCatlonal  contributions  to  the  transition  state.  Vardlaw  and  Marcus 

22 

and  Viswanathan  et  al.  have  investigated  these  issues  in  particular  for 
reactions  which  take  place  in  the  gas  phase.  Librational  motion, 
particularly  for  solution  phase,  is  important,  and,  indeed,  has  not  been 
neglected  entirely  in  our  work.  The  assumption  of  Bom-Oppenheimer 
separability  of  the  solute  vibrational  modes  from  degrees  of  freedom  of  the 
solvent  shows  how  specific  account  of  the  librations  is  included.  The  issue 
is  the  matter  of  time-scales.  The  Metropolis  samples  we  made  primarily 
involved  the  displacement  of  solvent.  The  solute  was  then  adjusted  to  the 
solvent  by  a  steepest  descent  calculation.  The  validity  of  this  sequence  of 
steps  rests  on  the  assumption  of  the  adiabatic  separation  of  the  solute 
intramolecular  vibrations  from  the  solvent  center-of-mass  motions;  we  did 
not  consider  the  intermolecular  coupling  of  internal  molecular  vibrations, 
although  it  is  possible  to  do  so.  We  assumed,  however,  that  the 
center-of-mass  motions  of  the  solute  were  similar  (in  scale)  from  those  of 
the  solvent.  Therefore,  displacement  of  the  solvent  according  to  the 
Metropolis  method  should  evolve  on  the  same  time-scale  as  libration  of  the 
solute.  Librational  contributions  are  then  automatically  included- -but  not 
distinguishable  from  solvent  contributions --in  the  averages. 

Librations  are  included  in  the  determination  of  the  tunnel  transition 
probability.  This  is  apparent  in  the  analysis  which  separates  the 
intramolecular  solute  vibrations  from  the  remaining  degrees  of  freedom  of 
the  solvent  and  center-of-mass  motion  of  the  solute.  This  is  implicit  (and 
can  be  made  explicit  with  a  vibrational  analysis  in  the  harmonic  limit)  in 
eq  (28) .  Eq  (45) ,  in  fact  must  be  summed  over  all  the  modes  which 
contribute.  For  a  system  of  limited  extent,  this  is  in  principle  possible. 
However,  even  for  the  size  of  system  which  we  have  considered  here,  the 
determination  of  the  eigenfrequencies  for  the  vibrations  together  with  the 


7» 

i 


sequence  of  steps  for  the  Metropolis  sampling  would  be  very  costly  in 
computer  time.  This  is  one  of  the  reasons  why  we  opted  to  use  the  simpler 
expression  eq  (47)  in  this  work. 

3t3 

Tapia  and  co-workers  have  recently  examined  an  actual  reactive  system 

in  solution:  the  protonation  of  3  methyl-but-l-yne-3-ol  in  water.  In  this 

work,  however,  they  did  not  optimize  the  transition  state  in  solution  as  we 

have  done.  Instead,  they  examined  the  role  of  solvent  in  the  system  in  its 

initial,  transition  and  final  states’,  the  transition  state  was  frozen  at  the 

vacuum  configuration.  They  indeed  found  that  there  were  important  solvent 

effects.  They  were  not  able,  however,  to  investigate  how  the  presence  of 

the  solvent  can  alter  the  structure  and  dynamics  of  the  solute  in  the 

transition  state  in  solution  over  that  determined  in  vacuum.  There  is 

growing  evidence2*  that  the  interaction  of  solute  with  solvent  and  with 

metal  surfaces  in  many  cases  alters  the  bonding  in  a  substrate  sufficiently 

to  see  in  the  vibrational  spectrum.  Although  a  study  of  type  of  Tapia  et 
23 

al.  can  yield  interesting  information  about  solvation  and  large  scale 
changes  in  solvation  in  going  from  the  initial  to  the  transition  state,  it 
is  probably  not  possible  to  obtain  accurate  estimates  of  other  quantities 
associated  with  the  transition  state.  Ve  think  that  our  approach  may  point 
to  way  to  filling  that  need. 

The  Metropolis  sampling  algorithm  is  not  the  only  possible  means  of 

modelling  chemical  reactivity.  Recently,  Wilson  and  Hynes  and 
23 

their  co-workers  reported  an  extensive  examination  of  the  type  of 
reaction  using  molecular  dynamics.  Their  approach  is  complimentary  to  ours 
in  the  sense  that  they  also  seek  to  follow  the  course  of  the  reaction  by 
employing  adequate  potential  energy  functions  for  the  interactions. 

Chemical  reactions  are  manifestly  time -dependent  phenomena,  and,  as  such,  it 
seems  immediately  clear  that  molecular  dynamics  is  the  most  natural  model  to 
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express  this  fact.  On  tha  other  hand,  the  transition  state  theory  of 
reactivity  aakes  extensive  use  of  ensemble  averages  which  can  be  examined  in 
an  essentially  time -independent  frame  of  reference.  In  this  sense,  the 
Metropolis/Monte  Carlo  may  prove  useful.  This  is  our  view. 

Pursuing  the  comparison  further,  Vilson  et  al.  used  their  molecular 
dynamics  algorithm  to  construct  the  potential  energy  surface  for  the 
reaction.  At  the  same  time,  however,  they  also  examined  the  frequency  of 
first  passage,  reflection  and  subsequent  passage.  Thus,  they  were  able  to 
estimate  the  transfer  frequency  from  two  points  of  view:  one  direct  and  the 
other  through  an  analysis  of  passage  over  and  tunnelling  through  the  average 
barrier.  There  are  issues  which  seem  better  handled  with  the  molecular 
dynamics  approach,  viz. ,  the  issue  of  "dynamic  friction"  and  the  question  of 
the  role  of  momentum  in  the  approach  to  and  passage  over  the  barrier.  In 
our  approach,  through  the  Monte  Carlo  algorithm,  we  seek  to  account  for  some 
similar  factors  through  the  detailed  analysis  of  the  tunnel  frequency.  It 
remains  to  be  settled  just  how  close  the  two  approaches  will  come  to  each 
other.  In  general,  the  Monte  Carlo  approach  is  computationally  cheaper;  it 
is  expected,  therefore,  that  one  day  it  may  be  useful  in  approximate,  but 
still  reasonably  accurate  and  predictive,  calculations  which  may  become 
generally  available. 

No  matter  how  simple  the  system,  the  examination  of  reactivity  consumes 
great  amounts  of  computer  time,  even  on  the  fastest  machines.  We  need  to 
look  for  more  direct,  and  probably  approximate,  methods  to  handle  the 
solvent.  As  an  alternative  to  scanning  for  meaningful  minima,  Warshel 
proposed  several  schemes,  one  of  which  is  his  "surface  constrained  soft 
sphere  dipole"  model  for  use  with  polar  solvents.  This  model  is  flexible 
and  accurate,  when  applied  to  the  determination  of  solvation  enthalpies.  As 
Varshel  notes,  the  enthalpy  is  the  leading  term- -the  entropy  is  a  second 
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order  effect. 
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It  would  seen  that  Uarshel's  approach  might  apply  to  the  determination 
of  transition  states  and  reaction  rates  as  we  have  proposed  here.  There  is, 
however,  a  problem  which  needs  to  be  resolved  before  Warshel's  model  of  the 
solvent  can  be  combined  with  our  approach.  Varshel's  technique  is  to  build 
shells  of  solvent  about  the  reactive  solute.  The  innermost  shell  is  the 
most  sensitive  to  states  and  changes  of  states  of  the  solute.  The  next 
layer  of  solvent  is  constrained  to  the  configuration  of  bulk  solution.  The 
states  of  the  solute  and  strongly  bound  solvent  can  then  be  determined  by 
means  of  steepest  descents  optimization,  and,  for  the  transition  state, 
saddle  point  optimizations.  The  difficulty  with  this  approach  is  the  fact 
that  static  optimization  of  this  type  only  finds  one  most  probable 
configuration.  Such  a  configuration  may  correspond  either  to  a 
configuration  which  is  completely  adiabatically  attained  or  to  one  for  which 
a  tunnel  transfer  is  required.  The  balance  between  tunnel  transfers  and 
over-the-top  passes  is  not  possible. 

Some  form  of  sampling  of  the  neighborhood  of  the  transition  state  is 
needed  in  order  to  determine  an  accurate  average  value  of  the 
pre- exponential  frequency  factor.  Warshel  has  already  noted  that  it  is 
possible  to  use  Metropolis/Monte  Carlo  methods  to  sample  states  in  the 
neighborhood  of  the  optimally  determined  state,  as  found,  for  example,  with 
the  use  of  a  steeptes  descent  routine.  In  this  way,  one  may  be  able  to 
determine  accurate  values  of  the  frequency  factor  with  a  modest  expense  in 
computer  time. 


Conclus ion 


We  have  suggested  that  the  reaction  rate  constant  can  be  constructed  as 


-33- 


an  ensemble  average  of  configurations  of  the  transition  state  which  are 
generated  according  to  a  modified  Metropolis  sampling  scheme .  For  each 
element  of  the  collection  of  configurations  in  the  ensemble  of  transition 
state  configurations,  the  transfer  species  either  occupies  a  local  energy 
well  or  must  tunnel  through  a  residual  barrier.  A  frequency  factor  can  be 
determined  for  each  case.  The  average  of  the  frequency  factors  associated 
with  these  two  possibilities  in  principle  yields  an  accurate  value  of  the 
pre- exponential  frequency  factor  in  the  Arrhenius  expression  for  the  rate 
constant.  Ve  believe  that  the  algorithm  which  we  suggest  shows  some 
potential  in  helping  us  to  examine  complicated  chemical  reactions  in 
solution. 

This  work  has  been  supported  in  part  by  contracts  with  the  Office  of 
Naval  Research,  Arlington,  VA. 

Appendix 

Basis  and  State  Functions  for  Use  in  Tunnelling  Calculations 

In  this  appendix,  we  outline  a  general  method  for  determining  state 
functions  which  can  be  used  in  tunnelling  calculations.  In  order  to 
construct  a  usable  Born-Oppenheimer  separation  of  the  local  modes  from  those 
of  the  surroundings,  it  is  necessary  to  develop  first  an  analysis  of  the 
reactive  subsystem.  Given  a  set  of  atomic  and  molecular  species  which  make 
up  the  reactive  molecular  system,  and  given  a  set  of  realistic  interaction 
potential  energy  functions,  the  first  objective  is  to  develop  an  analysis  of 
the  motions  of  the  isolated  molecular  system.  In  order  to  be  able  to 
follow  the  tunnelling  of  the  transfer  species,  it  is  necessary  to  identify 
the  parts  of  the  reactive  subsystems  and  then  to  connect  the  separate  parts 


in  an  analysis  which  focuses  on  the  terns  which  account  for  the  transfer. 

Vith  special  reference  to  the  transfer  species,  there  are  three 
principal  locations  for  this  species:  the  initial  configuration,  the 
location  of  the  transfer- species  in  the  transition  state,  and  finally  the 
location  of  the  species  in  the  final  state  configuration.  If  the  trajectory 
for  the  transfer- species  is  a  straight  line,  it  is  possible  to  consider  the 
tunnel  transition  only  with  reference  to  the  intial  and  final  stationary 
configurations.  If  the  reaction  path  is  geometrically  non-linear,  or 
oblique,  it  is  necessary  to  consider  a  piece-wise  linear  transfer;  precedent 
for  this  of  course  is  the  construction  of  basis  sets  of  atomic  orbitals  in 
the  familiar  LCAO  versions  of  the  molecular  orbital  theory.  An  illustration 
of  an  oblique  path  for  tunnelling  is  shown  in  Figure  4  for  a  hypothetical 
diffusion  of  a  species  from  one  vacancy  to  another. 

Although  the  initial  and  final  states  are  mechanically  stable,  the 
third  location,  the  transition  state,  is  not.  Nevertheless,  this  state  can 
be  included  in  a  linear  combination  in  order  to  connect  the  initial  with  the 
final  state  along  a  piece-wise  linear  path. 

We  imagine  a  dividing  plane  in  the  transition  state.  Left/initial  and 
right/final  states  are  defined  with  respect  to  this  plane.  We  next  consider 
a  set  of  harmonic  oscillator  basis  functions  which  are  centered  at  the 
initial,  transition,  and  final  states.  Let  the  set  be  (u^Cr^))  where  the 
index  i  is  a  function  label  and  the  index  I  is  a  label  for  the  point  of 
expansion.  These  functions  are  locally  orthonormal  with  respect  to  the 
point  of  expansion.  The  form  of  the  particular  functions  u^(r)  is  not 
important  at  this  stage,  but  they  are  typically  the  spherical  harmonic 
oscillator  or  similar  functions. 

We  construct  orthonormal  functions  for  the  initial  and  final  states  on 
either  side  of  the  dividing  plane.  These  functions  will  be  the  basis 


functions  for  the  construction  of  a  set  of  state  wavefunctions  for  the 


system.  The  elements  in  the  expansion  must  contain  functions  which  are 
centered  at  the  initial  (or  final)  state  and  the  transition  state;  the 
method  has  already  been  discussed  for  the  case  of  proton  transfer.  Thus, 
it  is  possible  to  allow  for  the  transitory  occupation  of  the  transition 
state  region  by  the  particle.  Functions  of  an  orthogonal  set  can  be  written 
as 


x  _ 

l(R)J  LCRJJ 


Id  u  (r  ) 
**  11  i  i 


(Al) 


where  the  coefficients  of  expansion  d are  determined  by  the  Gram- Schmidt 
method.  The  label  L(R)  in  eq  (Al)  indicates  a  function  for  the  left  (right) 
hand  side  of  the  transfer  hyperplane.  The  functions  6  are  now  used  to 
construct  the  state  functions  for  the  motion  of  the  particle  on  the  left 
(right)  hand  side  of  the  transfer  plane: 


^L(R)  “  ^  aL(R)j^L(R)j 

where  the  index  n  is  an  energy  level  index. 

Eigenvalues  and  eigenvectors  for  the  initial  and  final  transfer  states 
are  determined  on  the  left  and  right  hand  sides  of  the  transfer  plane;  this 
is  done,  as  usual,  by  means  of  a  variational  calculation.  For  H  the 
local  Hamiltonian  operator  for  the  motion  of  the  transfer  species  on  the 
left  (right)  of  the  transfer  plane,  a  set  of  eigenvalues  {«  .  }  results 

which  are  associated  with  the  migratory  species. 

As  long  as  the  basis  set  only  contains  functions  that  allow  for  motion 
of  the  transfer  species  between  two  stationary  points,  the  initial  and  the 
transition  state  or  the  transition  and  the  final  state,  these  functions 
alone  cannot  account  for  the  complete  migration. 

Next  the  optimally  determined  local  states  ^  are  used  to  construct 
state  functions  for  the  whole  transfer  system:  initial,  transition  and  final 


states.  Orthonormal  L/R  functions  are  constructed  by  Gram-Schmidt 


orthogonallzatlon: 
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*UR>  u  ;jloo  l(R> 

M 

with  the  coefficients  D  determined  by  the  orthonormalization.  The 

MUR)  J 


(A3) 


coefficients  (I-1,R)  in  the  expansion 

*  -  I  CiXl  (A4) 

are  determined  by  a  variational  calculation. 

The  result  of  all  of  these  steps  is  the  expression  of  the  matrix 
elements  for  the  transfer  in  terms  of  local  functions  u  (r^  centered  about 
each  of  the  three  primary  points  for  the  transfer:  Initial,  transition  and 
final  states. 

As  an  example,  consider  the  special  case  of  a  ground  state  transfer. 

As  basis  functions,  we  use  the  ground  state  spherical  oscillator  functions: 

uQ(r)  -  (a/w)3/*exp[-  j  ar2]  (A5) 

The  orthonormal  expansion  functions  for  the  left  and  right  hand  sides  are 
*t1  -  u  (r  )  (A6) 

Ll  0  a 

and 
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where  S>t  is  the  overlap  between  site  a  (on  the  left)  and  t  in  the 
transition  state;  there  are  similar  terms  for  6  and  A  These  functions 

*R1  rR2 

now  can  be  combined  to  get  the  left  and  right  hand  side  state  functions: 


4  -  a  4  +  a  4 

L  ll  Ll  L2L2 

4  -  a  4  +  4 

R  R1R1  R2R2 


(A8) 


The  lowest  eigenvalue  for  this  set  of  basis  functions  is  e 
construct  the  left  and  right  orthonormal  basis  functions: 


L(R) 


We  now 
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and  the  overlap  is  expressed  in  terms  of  the  basis  functions  u  (r^.  The 


system  energies  are  now  simply 


e  -  K.  •  v 1  !|(h*  -  v2  +  4v*'2 


(AlO) 


where  the  matrix  elements  have  an  obvious  origin.  These  matrix  elements 
also  can  be  expressed  in  terms  of  integrals  of  the  basis  set  u^r^) . 

We  now  have  a  representation  in  which  the  transfer  species  is  described 
in  an  initial  or  final  state  (left  or  right  state)  which  includes  an  account 
of  the  link  to  the  transition  state.  The  expansion  of  the  matrix  element 


in  eq  (AlO),  for  example,  indicates  a  connection  from  the  initial  to  the 


final  state  through  the  transition  state. 

There  has  been  some  discussion  of  this  type  of  approach  from  the  point 
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of  view  of  the  Bardeen  transfer  matrix  formalism.  The  approach  we  advocate 
here,  however,  is  slightly  different.  We  assume  that  there  exists  a  set  of 
well  behaved  and  denumerable  potential  energy  functions  which  operate 
between  the  species  in  the  system.  We  also  assume  that  the  structure  of  the 
system  is  known  or  at  least  specified.  The  transfer  formalism,  initiated 
by  Oppenheimer,30  works  well  for  barriers  which  are  not  fully  described, 
i.e.,  barriers  of  arbitrary  extent  and  height.  The  transfer  system,  as  we 
see  it,  is  as  well  characterized  as  is  the  atomic  substructure  in  a 
molecular  orbital  calculation.  We  simply  partition  the  wavefunctions  in 
such  a  manner  as  to  reveal  the  interactions  and  matrix  elements. 
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Table  Z 

Potential  energy  functions  and  parameters 


Bond 

Potential 

eA 

Characteristic  distance 

A-B 

Lennard-Jones 

e  - 

2.5 

a  -  1.41  A 

A-A 

Harmonic 

k  - 

8(80) 

r  -  2.75 

0 

A-S 

Lennard-Jones 

e  - 

0.05 

a  -  2.8 

B-S 

Lennard-Jones 

e  - 

0.05 

2  -  2.8 

The  form  of  the  Lennard-Jones  potential  used  is 
Vu  -  « (a/r)8[ (a/r)6  -  2] 

12 

and  the  pair  dissociation  energy  e  is  expressed  as  erg  x  10  The  harmonic 
potential  is 

Vr)  - 1 k  <r  - r/ 

and  k  is  expressed  as  dyn  cm  1  x  10  \ 


Table  II 


a.  Computed  quantities  for  the  ABA  inversion  with  k  -  8  x  lO^dyn  cm 


Temperature  (K) 

<E  > 

0 

<E*> 

<u> 

adiabatic ity 

4 

-9.429 

-9.594 

0.334 

0.0308 

50 

-9.413 

-9.437 

0.368 

0.2298 

100 

-9.157 

-9.233 

0.190 

0.7218 

200 

-8.782 

-8.774 

0.232 

0.5862 

300 

-8.588 

-8.587 

0.240 

0.4652 

b.  Computed  quantities  for  the 

ABA  inversion  with 

W 

1 

00 

X 

105dyn  cm  1 . 

4 

-9.484 

-9.557 

0.773 

0.00001 

50 

-9.400 

-9.311 

0.405 

0.0077 

100 

-9.262 

-9.052 

0.280 

0.0145 

200 

-8.790 

-8.686 

0.263 

0.0891 

300 


-8.568 


-8.504 


0.266 


0.0884 


Figure  Captions 


1.  The  initial  and  transition  states  for  the  inversion  of  S02  constrained 
to  two  dimensions. 

2.  The  energy  of  linear  S02  in  the  transition  state  for  the  inversion  as  a 
function  of  the  SO  bond  distance.  Only  a  symmetric  stretch  is  considered. 
The  energy  is  calculated  with  the  use  of  the  Carter,  et  al.  potential  (ref. 
7). 

3.  The  behavior  of  the  second  derivative  of  the  potential  energy  [En  - 
d  E/d z  ]  with  respect  to  the  displacement  of  the  S-atom  along  the  line 
normal  to  the  0-0  line.  The  curve  labelled  a  is  the  behavior  due  to  the 
Rydberg  pair  potential.  The  curve  b  is  the  behavior  of  the  many-body  part 
(which  accounts  for  angle  restoring  forces).  Curve  c  is  the  sum  of  a  and  b. 
The  region  of  curve  c  above  zero  indicates  the  region  of  stability  for  the 
motion  of  S  normal  to  the  0-0  line  in  the  transition  state. 

4.  The  course  of  diffusion  of  an  atom  from  one  location  to  a  vacancy 
along  an  oblique  path.  The  transition  state  does  not  lie  on  the  straight 
line  from  the  initial  to  the  final  site.  Account  of  the  transition  state 
can  be  incorporated  in  the  analysis  by  mixing  basis  functions  for  the 
diffusing  species  at  the  location  t-s  with  basis  functions  for  the  intial 


and  final  locations. 


Initial  State  Transition  state 
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